########################################################################
# PLS 206 — Week 3 live coding in class
# Grey Monroe, Fall 2026
########################################################################

########################################################################
# 1. Read the data and plot it first
########################################################################

# tomato_trial.csv: 120 greenhouse tomato plants (simulated), 20 of each
# variety x irrigation combination.
#   variety:    Roma, Cherry, Beefsteak
#   irrigation: full, deficit
#   calcium:    leaf calcium (mg/g)
#   yield_kg:   fruit yield per plant (kg)
#   ber:        blossom-end rot on the plant (1 = yes, 0 = no)
#   fruits:     number of fruits per plant

# install.packages(c("ggplot2", "emmeans", "car"))
tomato <- read.csv("https://greymonroe.github.io/PLS_206/data/tomato_trial.csv")

str(tomato)
table(tomato$variety, tomato$irrigation)

library(ggplot2)
ggplot(tomato, aes(variety, yield_kg, colour = irrigation)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.6)) +
  theme_classic()

# the same data as boxes with the points on top
ggplot(tomato, aes(variety, yield_kg, fill = irrigation)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.75)) +
  theme_classic()

# A 0/1 column (like ber) is stored as a number. To use it as a group
# (a fill colour, or a factor in lm()), wrap it in factor():
ggplot(tomato, aes(variety, yield_kg, fill = factor(ber))) +
  geom_boxplot() +
  theme_classic()

# The routine for every model below:
#   1. plot the raw data and save the plot (p_... <- ggplot(...))
#   2. fit the model
#   3. get its fitted values and 95% CIs for a grid of x values
#   4. add them to the saved plot as new layers (p_... + geom_...)

########################################################################
# 2. Two groups: a dummy variable
########################################################################

p_irr <- ggplot(tomato, aes(irrigation, yield_kg)) +
  geom_jitter(width = 0.1, alpha = 0.4) +
  labs(y = "Yield (kg per plant)") +
  theme_classic()
p_irr

# The intercept is the mean of the reference group (first in the
# alphabet: "deficit"); the coefficient is the other group minus it.

m_irr <- lm(yield_kg ~ irrigation, data = tomato)
round(coef(summary(m_irr)), 3)
tapply(tomato$yield_kg, tomato$irrigation, mean)

# the same test as a t-test with equal variances
t.test(yield_kg ~ irrigation, data = tomato, var.equal = TRUE)

# Make "full" the reference, so the coefficient reads as the cost of
# deficit irrigation:

tomato$irrigation <- factor(tomato$irrigation, levels = c("full", "deficit"))
m_irr <- lm(yield_kg ~ irrigation, data = tomato)
round(coef(summary(m_irr)), 3)

# fitted values with 95% CIs: predict(interval = "confidence") gives
# fit, lwr, upr for each row of newdata
fit_irr <- data.frame(irrigation = c("full", "deficit"))
fit_irr <- cbind(fit_irr, predict(m_irr, newdata = fit_irr, interval = "confidence"))
fit_irr

# the raw plot again, with the model on top
p_irr +
  geom_pointrange(data = fit_irr, aes(y = fit, ymin = lwr, ymax = upr), colour = "red")

########################################################################
# 3. Three groups: ANOVA
########################################################################

p_var <- ggplot(tomato, aes(variety, yield_kg)) +
  geom_jitter(width = 0.1, alpha = 0.4) +
  labs(y = "Yield (kg per plant)") +
  theme_classic()
p_var

# Each coefficient compares one variety with the reference (Beefsteak).
# anova() asks whether the varieties differ at all; TukeyHSD() compares
# every pair.

m_var <- lm(yield_kg ~ variety, data = tomato)
round(coef(summary(m_var)), 3)
anova(m_var)
TukeyHSD(aov(m_var))

fit_var <- data.frame(variety = c("Beefsteak", "Cherry", "Roma"))
fit_var <- cbind(fit_var, predict(m_var, newdata = fit_var, interval = "confidence"))
p_var +
  geom_pointrange(data = fit_var, aes(y = fit, ymin = lwr, ymax = upr), colour = "red")

# effect size: the share of variation explained by variety
ss <- anova(m_var)$"Sum Sq"
ss[1] / sum(ss)

########################################################################
# 4. Two factors and their interaction
########################################################################

# Does deficit irrigation cost every variety the same?

p_int <- ggplot(tomato, aes(irrigation, yield_kg, colour = variety)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.5), alpha = 0.35) +
  labs(y = "Yield (kg per plant)") +
  theme_classic()
p_int

m_int <- lm(yield_kg ~ variety * irrigation, data = tomato)
anova(m_int)                    # the variety:irrigation row first
round(coef(summary(m_int)), 3)

# With the interaction in, irrigationdeficit is the cost of deficit
# irrigation in Beefsteak only (the reference variety), and
# varietyRoma is Roma vs Beefsteak under full irrigation only.
# emmeans gives the cost in every variety directly:

library(emmeans)
emmeans(m_int, ~ irrigation | variety)
pairs(emmeans(m_int, ~ irrigation | variety))

# the six fitted means with their CIs; lines that are not parallel are
# the interaction
fit_int <- expand.grid(variety = c("Beefsteak", "Cherry", "Roma"),
                       irrigation = c("full", "deficit"))
fit_int <- cbind(fit_int, predict(m_int, newdata = fit_int, interval = "confidence"))
p_int +
  geom_pointrange(data = fit_int, aes(y = fit, ymin = lwr, ymax = upr),
                  position = position_dodge(width = 0.5)) +
  geom_line(data = fit_int, aes(y = fit, group = variety),
            position = position_dodge(width = 0.5))

########################################################################
# 5. summary() vs anova()
########################################################################

# Same model, a different reference variety. The coefficients change;
# the anova() table does not.

tomato$variety2 <- factor(tomato$variety, levels = c("Roma", "Cherry", "Beefsteak"))
m_int2 <- lm(yield_kg ~ variety2 * irrigation, data = tomato)
round(coef(summary(m_int2)), 3)
anova(m_int2)
# irrigationdeficit is now the cost of deficit in Roma (0.74 kg), and the
# interaction rows compare Cherry and Beefsteak with Roma: Cherry loses
# about as much as Roma (p = 0.24), Beefsteak loses much more.
# The fitted values are identical to m_int's, so the plot is the same:
all.equal(fitted(m_int), fitted(m_int2))

########################################################################
# 6. Groups and a slope: ANCOVA
########################################################################

# Leaf calcium varies from plant to plant. One slope, a shift for
# deficit irrigation, then a separate slope for each irrigation level.

p_anc <- ggplot(tomato, aes(calcium, yield_kg, colour = irrigation)) +
  geom_point(alpha = 0.5) +
  labs(x = "Leaf calcium (mg/g)", y = "Yield (kg per plant)") +
  theme_classic()
p_anc

m_anc <- lm(yield_kg ~ calcium + irrigation, data = tomato)
round(coef(summary(m_anc)), 4)
# The deficit effect shrinks from -1.32 (section 2) to -0.93 kg here:
# deficit plants have less leaf calcium (about 3 mg/g less), so part of
# the irrigation effect was really a calcium effect.
tapply(tomato$calcium, tomato$irrigation, mean)

# fitted lines with 95% CI bands
fit_anc <- expand.grid(calcium = seq(7, 31, by = 0.5), irrigation = c("full", "deficit"))
fit_anc <- cbind(fit_anc, predict(m_anc, newdata = fit_anc, interval = "confidence"))
p_anc +
  geom_ribbon(data = fit_anc, aes(y = fit, ymin = lwr, ymax = upr, fill = irrigation),
              alpha = 0.2, colour = NA) +
  geom_line(data = fit_anc, aes(y = fit), linewidth = 1)

m_anc2 <- lm(yield_kg ~ calcium * irrigation, data = tomato)
round(coef(summary(m_anc2)), 4)
# calcium:irrigationdeficit is far from significant: the slopes are
# parallel, so keep the simpler model.
# Don't read the irrigationdeficit row of this model (p = 0.50): with
# the interaction in, it is the deficit effect at calcium = 0, and no
# plant has less than 7 mg/g (Monday's centring slides).

# the same raw plot with the separate-slopes model: nearly parallel
fit_anc2 <- expand.grid(calcium = seq(7, 31, by = 0.5), irrigation = c("full", "deficit"))
fit_anc2 <- cbind(fit_anc2, predict(m_anc2, newdata = fit_anc2, interval = "confidence"))
p_anc +
  geom_ribbon(data = fit_anc2, aes(y = fit, ymin = lwr, ymax = upr, fill = irrigation),
              alpha = 0.2, colour = NA) +
  geom_line(data = fit_anc2, aes(y = fit), linewidth = 1)

########################################################################
# 7. Unbalanced data: ANOVA types
########################################################################

# Drop 25 plants at random; the cells are no longer equal, so the order
# of terms in anova() matters. car::Anova(type = 2) adjusts each main
# effect for the other.

set.seed(7)
unbal <- tomato[-sample(nrow(tomato), 25), ]
table(unbal$variety, unbal$irrigation)

p_unbal <- ggplot(unbal, aes(irrigation, yield_kg, colour = variety)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.5), alpha = 0.35) +
  labs(y = "Yield (kg per plant)") +
  theme_classic()
p_unbal

anova(lm(yield_kg ~ variety + irrigation, data = unbal))
anova(lm(yield_kg ~ irrigation + variety, data = unbal))
library(car)
m_unbal <- lm(yield_kg ~ variety + irrigation, data = unbal)
Anova(m_unbal, type = 2)

# The order of terms changes the anova() table, not the model: the
# fitted values (and so the plot) are the same either way.
fit_unbal <- expand.grid(variety = c("Beefsteak", "Cherry", "Roma"),
                         irrigation = c("full", "deficit"))
fit_unbal <- cbind(fit_unbal, predict(m_unbal, newdata = fit_unbal, interval = "confidence"))
p_unbal +
  geom_pointrange(data = fit_unbal, aes(y = fit, ymin = lwr, ymax = upr),
                  position = position_dodge(width = 0.5)) +
  geom_line(data = fit_unbal, aes(y = fit, group = variety),
            position = position_dodge(width = 0.5))

########################################################################
# 8. A yes/no response: blossom-end rot
########################################################################

# Blossom-end rot is a calcium disorder. The coefficients are on the
# log-odds scale; plogis() turns log-odds into a probability, and exp()
# turns a coefficient into an odds ratio.

table(tomato$ber, tomato$irrigation)

p_ber_irr <- ggplot(tomato, aes(irrigation, ber)) +
  geom_jitter(width = 0.1, height = 0.04, alpha = 0.4) +
  labs(y = "Blossom-end rot (0 = no, 1 = yes)") +
  theme_classic()
p_ber_irr

m_ber_irr <- glm(ber ~ irrigation, family = binomial, data = tomato)
round(coef(summary(m_ber_irr)), 3)
plogis(coef(m_ber_irr)[1])           # P(rot), full irrigation
plogis(sum(coef(m_ber_irr)))         # P(rot), deficit
prop.table(table(tomato$ber, tomato$irrigation), 2)
exp(coef(m_ber_irr))                 # odds ratio for deficit
exp(confint(m_ber_irr))

# Fitted values and CIs for a GLM: build the interval on the link scale
# (log-odds here), where it is symmetric, then turn all three numbers
# back into the response scale (probability) with the inverse link.
# Works for any glm(); used for every GLM below.
glm_ci <- function(model, newdata) {
  p <- predict(model, newdata = newdata, type = "link", se.fit = TRUE)
  inv <- family(model)$linkinv          # plogis for binomial, exp for poisson
  newdata$fit <- inv(p$fit)
  newdata$lwr <- inv(p$fit - 1.96 * p$se.fit)
  newdata$upr <- inv(p$fit + 1.96 * p$se.fit)
  newdata
}
fit_ber_irr <- glm_ci(m_ber_irr, data.frame(irrigation = c("full", "deficit")))
fit_ber_irr

# red: the fitted probability of rot with its 95% CI
p_ber_irr +
  geom_pointrange(data = fit_ber_irr, aes(y = fit, ymin = lwr, ymax = upr), colour = "red")

# now calcium as the predictor
p_ber <- ggplot(tomato, aes(calcium, ber)) +
  geom_jitter(width = 0, height = 0.04, alpha = 0.4) +
  labs(x = "Leaf calcium (mg/g)", y = "Blossom-end rot (0 = no, 1 = yes)") +
  theme_classic()
p_ber

m_ber <- glm(ber ~ calcium, family = binomial, data = tomato)
round(coef(summary(m_ber)), 3)
exp(coef(m_ber)["calcium"])          # odds ratio per mg/g of calcium

# predicted probability of rot at three calcium levels
new <- data.frame(calcium = c(10, 20, 30))
predict(m_ber, newdata = new, type = "response")

# red: the fitted probability curve with its 95% CI band
fit_ber <- glm_ci(m_ber, data.frame(calcium = seq(7, 31, by = 0.5)))
p_ber +
  geom_ribbon(data = fit_ber, aes(y = fit, ymin = lwr, ymax = upr), fill = "red", alpha = 0.2) +
  geom_line(data = fit_ber, aes(y = fit), colour = "red", linewidth = 1)

########################################################################
# 9. How often is it right?
########################################################################

# Compare the accuracy with always guessing the most common answer.

pred <- as.numeric(fitted(m_ber) > 0.5)
table(predicted = pred, observed = tomato$ber)
mean(pred == tomato$ber)
mean(tomato$ber == 0)

# Calcium alone gets about 72% right, against 56% for always guessing
# "no rot". A clear effect (look at its p-value) still misclassifies
# more than a quarter of the plants: a real effect and a good
# classifier are different things.

########################################################################
# 10. Deviance, then several predictors
########################################################################

# What is deviance? For each plant, take the probability the model gave
# to what actually happened (p if it rotted, 1 - p if it didn't). A good
# model gives a high probability to what happened. Each plant's penalty
# is -2 * log(that probability): near 0 if the model was confident and
# right, large if it was confident and wrong. Deviance is the sum of the
# penalties, the GLM's version of the residual sum of squares.

tomato$p_fit <- fitted(m_ber)
tomato$p_happened <- ifelse(tomato$ber == 1, tomato$p_fit, 1 - tomato$p_fit)
tomato$penalty <- -2 * log(tomato$p_happened)
sum(tomato$penalty)
deviance(m_ber)                      # the same number

# The null deviance is the same sum for a model with no predictor: every
# plant gets the overall rot rate.
p0 <- mean(tomato$ber)
sum(-2 * log(ifelse(tomato$ber == 1, p0, 1 - p0)))
m_ber$null.deviance                  # the same number

# The picture: a line runs from each plant to the model's curve, down
# from 1 for plants that rotted, up from 0 for plants that didn't.
# Short line = the model gave a high probability to what happened.
# (No jitter here, so each line starts exactly at 0 or 1.)
p_dev <- ggplot(tomato, aes(calcium, ber)) +
  geom_point(alpha = 0.4) +
  labs(x = "Leaf calcium (mg/g)", y = "Blossom-end rot (0 = no, 1 = yes)") +
  theme_classic()
p_dev +
  geom_segment(data = tomato, aes(xend = calcium, y = ber, yend = p_fit, colour = factor(ber)), alpha = 0.6) +
  geom_line(data = fit_ber, aes(y = fit), colour = "red", linewidth = 1) +
  labs(colour = "rot", title = paste("model: deviance =", round(deviance(m_ber), 1)))
p_dev +
  geom_segment(data = tomato, aes(xend = calcium, y = ber, yend = p0, colour = factor(ber)), alpha = 0.6) +
  geom_hline(yintercept = p0, colour = "red", linewidth = 1) +
  labs(colour = "rot", title = paste("no predictor: null deviance =", round(m_ber$null.deviance, 1)))

# each plant's penalty: the confidently wrong plants dominate
ggplot(tomato, aes(calcium, penalty, colour = factor(ber))) +
  geom_point() +
  labs(x = "Leaf calcium (mg/g)", y = "Penalty, -2 log(probability of what happened)",
       colour = "rot") +
  theme_classic()

# Calcium and variety together. (Irrigation acts mostly through calcium,
# see section 6, so we leave it out here.)

p_ber2 <- ggplot(tomato, aes(calcium, ber, colour = variety)) +
  geom_jitter(width = 0, height = 0.04, alpha = 0.5) +
  labs(x = "Leaf calcium (mg/g)", y = "Blossom-end rot (0 = no, 1 = yes)") +
  theme_classic()
p_ber2

m_ber2 <- glm(ber ~ calcium + variety, family = binomial, data = tomato)
round(coef(summary(m_ber2)), 3)
round(exp(coef(m_ber2))[-1], 2)                 # odds ratios (skip the intercept)
1 - deviance(m_ber2) / m_ber2$null.deviance     # pseudo-R2
anova(m_ber, m_ber2, test = "Chisq")            # does variety add anything?

# one fitted curve per variety, same slope on the log-odds scale
fit_ber2 <- glm_ci(m_ber2, expand.grid(calcium = seq(7, 31, by = 0.5),
                                       variety = c("Beefsteak", "Cherry", "Roma")))
p_ber2 +
  geom_ribbon(data = fit_ber2, aes(y = fit, ymin = lwr, ymax = upr, fill = variety),
              alpha = 0.15, colour = NA) +
  geom_line(data = fit_ber2, aes(y = fit), linewidth = 1)

# accuracy of the bigger model, same code as section 9
pred2 <- as.numeric(fitted(m_ber2) > 0.5)
mean(pred2 == tomato$ber)

########################################################################
# 11. A count: fruits per plant
########################################################################

p_fruit <- ggplot(tomato, aes(variety, fruits, colour = irrigation)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.5), alpha = 0.35) +
  labs(y = "Fruits per plant") +
  theme_classic()
p_fruit

# Poisson regression: the coefficients are on the log scale; exp() turns
# them into rate ratios (how many times more fruit).

m_fruit <- glm(fruits ~ variety + irrigation, family = poisson, data = tomato)
round(coef(summary(m_fruit)), 3)
round(exp(coef(m_fruit)), 2)

# fitted mean counts with 95% CIs (glm_ci uses exp() for the log link)
fit_fruit <- glm_ci(m_fruit, expand.grid(variety = c("Beefsteak", "Cherry", "Roma"),
                                         irrigation = c("full", "deficit")))
p_fruit +
  geom_pointrange(data = fit_fruit, aes(y = fit, ymin = lwr, ymax = upr),
                  position = position_dodge(width = 0.5))

# Overdispersion check: about 1 is fine, well above 1.5 is a problem
# (the lecture's aphid counts gave 5.7 and needed MASS::glm.nb()).

deviance(m_fruit) / df.residual(m_fruit)

########################################################################
# 12. Report it
########################################################################

# An effect size with its CI, on a scale a reader understands.

exp(coef(m_ber)["calcium"])
exp(confint(m_ber)["calcium", ])
exp(coef(m_fruit)["irrigationdeficit"])
exp(confint(m_fruit)["irrigationdeficit", ])